% This file (1) creates a map of WB investments across LAC from 2005-2020
% and (2) maps WB investments to country grid cells for a comparison of 
% model-based investments to WB investments
% common parameters to all graphs
bright_line = 0.5;  % higher=brighter
bright_nodes = 0.4; % higher=brighter
adj_width = 6;
max_bright = 1.5;
markersize = 5;
color_growth = [ 0 1 0 ];   % red = [ 1 0 0 ]  green = [ 0 1 0 ] blue = [ 0 0 1 ]
color_shrink = [ 1 0 0 ];

addpath('../../Toolbox/export_fig/');

folders;

country_list_short;

% get WB projects -- total amount by coordinates
wbproj = shaperead(['../../data/generated/wb_api_bycoords.shp']);

% For each location, the funding is a constant share of total/number of places
country_borders = shaperead([path_raw_mapdata, '/COUNTRY_SHP/countries_shp/countries.shp']);

% filter out shapefile to rows that are the countries in our analysis
keep = zeros(length(country_borders), 1);
for j = 1:length(country_borders)
    if ismember(country_borders(j).ISO2, country_names(:, 2))
        keep(j) = 1;
    end
end

% Plot remaining countries 
country_borders = country_borders(keep==1);

% Draw plot
amount = cell2mat({wbproj.amount_y}');

mapshow(country_borders, 'FaceColor', 'white', 'EdgeColor', 'black')
hold on;
for i = 1:length(wbproj)
    if amount(i) == 0
        mapshow(wbproj(i), 'Marker','o', ...
            'MarkerFaceColor', 'red', 'MarkerSize', 3);
    else
         mapshow(wbproj(i), 'Marker','o', ...
            'MarkerFaceColor', 'none', 'MarkerEdgeColor', 'red', 'MarkerSize', ...
            3+amount(i)/max(amount)*25);
    end
    hold on;
end

margin = 0.5;
Xlim = [ min( cell2mat({country_borders.X}) )-margin;max( cell2mat({country_borders.X}) )+margin ];
Ylim = [ min( cell2mat({country_borders.Y}) )-margin;max( cell2mat({country_borders.Y}) )+margin ];
set( gca,'XTick',[],'YTick',[],'XLim',Xlim,'YLim',Ylim,'Box','on' )

export_fig([ path_notes_images,  'wb_projects.eps' ], '-depsc');

% load in the grid cells for each country and aggregate the (1) number of
% projects and (2) amount of the projects in each cell 
% compare with the investment in each grid cell per country

corrs_count = zeros(Ncountries, 1);
corrs_amount = zeros(Ncountries, 1);
corrs_population = zeros(Ncountries, 1);
corrs_population2 = zeros(Ncountries, 1);

list_of_correlations = table(country_names(:, 3), corrs_amount, corrs_count, corrs_population ,corrs_population2);
for country_n=1:Ncountries

    % load in expansion file name
    country_icc = char( countries(country_n) );  % country ICC code
    country = icc2name( country_icc );
    
    [ country ]
    x_ver_hor = get_threshold(country);
    x_diag = x_ver_hor;
    x_thresh =  x_ver_hor;
    
    default_cellsize = 0.5;

    % load in the counterfactual grid 
    filename = [path_save_cfactuals, country, '_diag', num2str( x_diag ), '_hor', num2str( x_ver_hor ), ...
                '_a1_rho0_alpha0.4_sigma5_beta0.13_gamma0.1_nu1_mobil0_cong1_ngoods10_cellsize0.5_WP0_osm0_cfactual_exp_eng.mat'];
    load(filename)
    
    % load in the actual grod 
    load( [ path_save_grids,country,'_grid_',...
        num2str( x_diag ),'_',...
        num2str( x_ver_hor ), '_',...
        num2str( default_cellsize ),'_connected1_WP0_osm0.mat' ] );

    % grab the parts of this we need 
    gridmap = country_graph.gridmap;
    g = cfactual.g;
    results_actual = cfactual.results_actual;
    results_exp = cfactual.results_cfactual;
    unique_edges=country_graph.unique_edges;
    country_bounds=country_graph.country_bounds;
    discretized_roads=country_graph.discretized_roads;
    places2=country_graph.places2;
    
    sum_avI_obs = sum( g.avI,2 );
    sum_avI_exp = sum( results_exp.Ijk,2 );
    
    % lanes in calibrated and counterfactual model
    I_obs = zeros( length(unique_edges),1 );
    I_exp = zeros( length(unique_edges),1 );
    for j=1:length(unique_edges)
        I_obs(j) = g.avI( unique_edges(j,1),unique_edges(j,2) );
        I_exp(j) = results_exp.Ijk( unique_edges(j,1),unique_edges(j,2) );
    end

    dI = I_exp-I_obs;
    dI(dI>prctile(dI, 90)) = prctile(dI, 90);

    % Compute total investment at the gridcell level and get the
    % corresponding level of $$ Wb money spent

    for j=1:length(gridmap)
        gridmap(j).wbproj_amount = 0;
        gridmap(j).wbproj_count = 0;
        gridmap(j).dI = sum_avI_exp(j)-sum_avI_obs(j);
        gridmap(j).proj_num = [];
        
        % check which projects (+ how many) are inside the gridcell
        for k=1:length(wbproj)
            if inpolygon(wbproj(k).Y, wbproj(k).X, gridmap(j).Y, gridmap(j).X)
                gridmap(j).wbproj_amount = gridmap(j).wbproj_amount + wbproj(k).amount_y;
                gridmap(j).wbproj_count =  gridmap(j).wbproj_count + 1;   
                gridmap(j).proj_num = [ gridmap(j).proj_num k];
            end
        end        
    end
    
    wb_count = cell2mat({gridmap.wbproj_count})';
    wb_amount = cell2mat({gridmap.wbproj_amount})';
    model_investment = cell2mat({gridmap.dI})';
    population = cell2mat({gridmap.population})';
    
    % get correlations across 
    rho_count = corr(model_investment,  wb_count);
    rho_amount = corr(model_investment,  wb_amount);
    rho_pop = corr(population,  wb_amount);
    rho_pop2 = corr(population,  wb_count);

    if sum(cell2mat({gridmap(:).wbproj_count }))>1
        list_of_correlations(strcmp(list_of_correlations.Var1, country) , 2) = num2cell(rho_amount);
        list_of_correlations(strcmp(list_of_correlations.Var1, country) , 3) = num2cell(rho_count);    
        list_of_correlations(strcmp(list_of_correlations.Var1, country) , 4) = num2cell(rho_pop);  
        list_of_correlations(strcmp(list_of_correlations.Var1, country) , 5) = num2cell(rho_pop2);    
    end
   
end

% write to csv 
writetable(list_of_correlations, ['../../../notes/wb_Iexp.csv', ])


            